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Abstract 

We study the lattice spacing dependence, or scaling, of physical quantities using the highly 
improved staggered quark (HISQ) action introduced by the HPQCD/UKQCD collaboration, com- 
paring our results to similar simulations with the asqtad fermion action. Results are based on 
calculations with lattice spacings approximately 0.15, 0.12 and 0.09 fm, using four flavors of dy- 
namical HISQ quarks. The strange and charm quark masses are near their physical values, and 
the light-quark mass is set to 0.2 times the strange-quark mass. We look at the lattice spacing 
dependence of hadron masses, pseudoscalar meson decay constants, and the topological suscepti- 
bility. In addition to the commonly used determination of the lattice spacing through the static 
quark potential, we examine a determination proposed by the HPQCD collaboration that uses the 
decay constant of a fictitious "unmixed ss" pseudoscalar meson. We find that the lattice artifacts 
in the HISQ simulations are much smaller than those in the asqtad simulations at the same lattice 
spacings and quark masses. 

PACS numbers: 12.38.Gc,14.20.Dh 
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I. INTRODUCTION AND MOTIVATION 



The "highly improved staggered quark", or HISQ, action was developed by the 
HPQCD/UKQCD collaboration to reduce the lattice artifacts associated with staggered 
quarks in lattice QCD calculations [IH3]. While significantly more expensive than the asq- 
tad action used in the MILC collaboration's long-running program of QCD simulations with 
three dynamical quark flavors it is still very economical compared with non-staggered 
quark actions. 

The initial studies of the HISQ action by the HPQCD/UKQCD collaboration demon- 
strated the reduction of taste symmetry breaking and improvements in the dispersion rela- 
tion for the charm quark by using the HISQ action for valence quarks on quenched lattices 
and lattices generated with asqtad sea quarks [TH3]. Further work with this action, again 
implemented for the valence quarks with asqtad sea quarks, has demonstrated impressive 
precision for charmonium and heavy-light meson physics [5H7]. 

As a first stage in a complete program of QCD simulations using the HISQ action for 
dynamical quarks, we have generated ensembles of lattices at three different lattice spacings 
with four flavors of dynamical quarks, where the light-quark mass is fixed at two-tenths of the 
strange quark mass, and the strange and charm quark masses are near their physical values. 
This allows us to test scaling, or dependence of calculated quantities on the lattice spacing. 
The purpose of this paper is to report on these tests at fixed quark mass. Where possible, 
we compare the lattice spacing dependence of physical quantities with the HISQ action to 
their dependence using the asqtad action at the same quark mass and lattice spacings. We 
look at the static quark potential, splittings among the different tastes of pions, masses of 
the rho and nucleon, pseudoscalar meson decay constants and the topological susceptibility. 
We emphasize that all of this is done at a fixed, and unphysically large, light-quark mass - 
our purpose here is to make a controlled study of the dependence on lattice spacing. 

II. METHODS AND LATTICE DATA 

There are four major differences between these HISQ simulations and our earlier asqtad 
simulations. 

First, the HISQ simulations include the effects of a dynamical charm quark. We expect 
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that the effects of dynamical charm will be very small for the quantities studied here, but 
with modern algorithms it is cheap to include the charm quark, and we plan to investigate 
quantities involving dynamical charm in the future. 

Second, the one-quark-loop contributions to the perturbative calculation of the coeffi- 
cients in the Symanzik improved gauge action are included. At the time the asqtad sim- 
ulation program was started these corrections were not available, but they have now been 
computed for both the asqtad and HISQ actions, and are unexpectedly large [8]. 

Third, in the HISQ action the parallel transport of quark fields is done with a link that 
is highly smeared. Specifically, it is first smeared using a "fat7" smearing, then projected 
onto a unitary matrix, and then smeared again with an "asqtad" smearing j3]. The use of 
the asqtad smearing in the second iteration, together with the addition of the Naik term, or 
third-nearest-neighbor coupling, in p, insures that the fermion action is formally order a 2 
improved. The use of two levels of smearing produces a smooth gauge field as seen by the 
quarks, and this explains the reduced taste symmetry violations. 

Finally, the third-nearest-neighbor term in the charm quark p is modified to improve 
the charm quark dispersion relation [3]. These last two differences combine to make up 
what is usually meant by "the HISQ action", although in principle they could be introduced 
independently. 

Where practical, since our purpose is to compare the lattice artifacts in the two actions, 
we use the same analysis for the HISQ data as was used for the asqtad data. 

Table [I] shows the parameters of the three HISQ runs used in these tests. Detailed 
information about the asqtad ensembles can be found in Ref. 

The HISQ lattices were generated using the rational hybrid Monte Carlo (RHMC) algo- 
rithm [9]. Issues with implementing this algorithm for the HISQ action have been discussed 
in Ref. [10J. We used different molecular dynamics step sizes for the gauge and fermion 
parts of the action, with three gauge steps for each fermion step [11]. We used the Omelyan 
integration algorithm in both the gauge and fermion parts [TTI IT2] . Five pseudofermion fields 
were used, each with a rational function approximation for the fractional powers. The first 
implements the ratio of the roots of the determinants for the light and strange sea quarks 
to the determinant for three heavy "regulator" quarks with mass am r = 0.2. That is, it 
corresponds to the weig ht det (M(m,)) 1/2 det (M(m s )) 1/4 det (M(m r ))" 3/4 . The next three 
pseudofermion fields each implement the force from one flavor of the regulator quark, or the 
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TABLE I: Parameters of the HISQ runs with mi = 0.2 m s . Here e^v is the correction for the three 
link (Naik) term in the charm quark action. These values differ slightly from the expression in 
Appendix [A] because they do not include the distinction between bare and tree-level quark mass 
(see Eq. (24) in Ref. [3].) The expression in Appendix |A| is used in all more recent ensembles. The 
number of equilibrated lattices is N[ ats . The separation of the lattices in simulation time is St, 
the length of a trajectory in simulation time is Lt, the molecular dynamics step size is e, and the 
fraction of trajectories accepted is "ace". Our definition of the step size is such that there is one 
evaluation of the fermion force per step, so a complete cycle of the Omelyan integration algorithm 
includes two fermion- action steps and six gauge-action steps. The physical lattice spacing given in 
this table uses the three flavor determination of n = 0.3117(6)(l3^) fm made using f n to set the 
scale on the asqtad ensembles [15]. It should be noted that when chiral and continuum limits of 
2+1+1 flavor calculations are completed, a 2+1+1 flavor determination of r\ will supercede this. 
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ami am s am c 


size 


M 


Ni a t s S t L t e acc. 


ri/a a (fm) 


5.8 
6.0 
6.3 


0.013 0.065 0.838 -0.3582 
0.0102 0.0509 0.635 -0.2308 
0.0074 0.037 0.440 -0.1205 


16 3 x 48 
24 3 x 64 
32 3 x 96 


0.85535 
0.86372 
0.874164 


1021 5 1.0 0.033 0.73 
1040 5 1.0 0.036 0.66 
878 6 1.5 0.031 0.68 


2.041(10) 0.1527(j^ 6 ) 
2.574(5) 0.1211(^ 2 ) 
3.520(7) 0.0886(^) 



fourth root of the corresponding determinant [13]. The final pseudofermion field implements 
the dynamical charm quark. 

Rational function approximations were used for the fractional powers of the matrices 
[H]. In the molecular dynamics evolution we used a 9'th order approximation for the 
pseudofermion field containing the light quarks, and a 7'tli order approximation for the 
three regulator fields and the charm quark pseudofermion. For the heat bath updating of 
the pseudofermion fields and for computing the action at the beginning and end of the 
molecular dynamics trajectory we used ll'th order and 9'th order approximations. These 
approximations comfortably exceeded the required accuracy, but since a multimass conju- 
gate gradient routine is used for the sparse matrix solutions, adding extra terms in these 
approximations has minimal cost. 

In order to make this paper self-contained, we summarize the action in Appendix |X| and 
discuss some algorithmic issues specific to the HISQ action in Appendices |B| |C| and |D} 
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FIG. 1: Autocorrelation C^t m simulation time of the plaquette (left panel), strange quark ipif) 
(center panel) and topological charge (right panel). Note that the horizontal scale is different in 
each of the three panels. Errors on the autocorrelation were estimated by dividing the time series 
into five subsets and averaging the autocorrelations from each subset. The vertical arrows in the 
left panel indicate the time separation between stored lattices, used in computing the potential, 
spectrum and other quantities. 

III. AUTOCORRELATIONS IN SIMULATION TIME 

Estimating statistical errors on any physical quantity requires taking into account the 
fact that successive sample configurations are not completely statistically independent, and 
calculations of statistical errors that ignore these autocorrelations are generally underesti- 
mates of the true errors. The amount of autocorrelation depends strongly on the quantity 
under consideration, so we present autocorrelations for a few simple but relevant quantities. 
To parameterize autocorrelations we use the dimensionless coefficient 



where Xi is the measurement at simulation time i and At is the time separation of the two 
measurements. As discussed above, in these simulations successive lattices were saved at 
time separations At = 5 for the a = 0.15 and 0.12 fm ensembles, and At = 6 for the 0.09 
fm ensemble. However, measurements of the plaquette and ipip were made every trajectory. 
Note that determination of these autocorrelation coefficients is numerically difficult, even 
on time series of order 1000 lattices. This is partly because of the practical necessity of 
using the average ((£»)) from our simulation, rather than the true average. (Note, however, 
that for the topological charge we know that the true average is zero.) Estimation of errors 
on these coefficients is also noisy. Here we have estimated the errors from the variance of 
autocorrelations measured on five separate segments of the time series, but for the central 
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TABLE II: Autocorrelation C&t of various quantities between successive lattices in the ensembles. 
Lattices are separated by five time units for a = 0.15 and 0.12 fm, and by six time units for 
a = 0.09 fm. As discussed in the text, the autocorrelations for tpip are between estimates made 
with one random source. Autocorrelations for the correlators {tt(Q)tt{D)) and (p(0)p(D)) are given 
at a spatial distance D which is the minimum distance used in a typical fit for the mass. For the 
pion correlator these distances are D = 15, 20 and 30 respectively, and for the p correlator they 
are D = 6, 7 and 10 respectively. For the pion and rho mass and the autocorrelations are from 
single elimination jackknife samples. 



operator 


0.15 fm 


0.12 fm 


0.09 fm 


□ 


0.311(25) 


0.300(10) 


0.359(14) 


flight 


0.135(24) 


0.151(27) 


0.192(34) 


^/'V'strange 


0.234(38) 


0.265(27) 


0.259(19) 


<7r(0)7r(£>)) 


0.034(40) 


0.084(46) 


0.177(21) 


(P(0)P(D)) 


0.055(24) 


0.074(24) 


0.061(18) 


m n 


0.008(14) 


0.182(35) 


0.249(51) 


U 


0.123(21) 


0.150(23) 


0.184(45) 


m p 


0.036(38) 


0.045(09) 


0.002(24) 


Qtopo 


na 


0.500(25) 


0.754(36) 



value quote the result from the full time series. 

The autocorrelations can be taken into account either by blocking the data (averaging over 
intervals of time) and then computing the average of the blocked values, or by multiplying 
the error estimate ignoring autocorrelations by the factor 

Ji+zjyt , (2) 

with the sum suitably truncated. For complicated functions of observables a jackknife anal- 
ysis can be used, and Eqs. ([T| and ^ applied to the sequence of jackknife results. 

We begin with the plaquette and strange quark x/jip, simple observables which were mea- 
sured at each trajectory. The first two panels of Figure [T] show the autocorrelation of these 
quantities as a function of separation in simulation time. Here ipip is estimated using a single 
random source vector. Thus part of its variance comes from the random source, and part 
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from the variation of the lattice. For this reason, its autocorrelation does not approach one 
at small time. We show the strange quark ifjif) because it generally shows larger autocor- 
relations than the light quark ipip. Also, relevant to future ensembles at other light quark 
masses, it will be useful to compare autocorrelations using ifjip at a fixed physical quark 
mass. These two simple quantities provide a good illustration of how autocorrelations differ 
among various quantities. 

Table |TT] shows these quantities at the time separation of the stored lattices, and a se- 
lection of autocorrelations of more physically relevant quantities. In particular, it contains 
autocorrelations of the pion and rho correlators ((ir(0)ir(D)) and (p(0) p(D))) at a distance 
D equal to the minimum distance that might be used in a mass fit, and would be one of 
the important contributors to the mass. This table also contains autocorrelations of single 
elimination jackknife measurements of the pion mass, the pion decay constant (amplitude 
of a pion correlator), and the rho meson mass. 

The topological charge is generally expected to have a long autocorrelation time. In 
fact, in the continuum limit tunnelings would be expected to be completely suppressed in 
a simulation algorithm where the configurations evolve continuously. Such a simulation 
would still give correct results in infinite volume, but would have power-law finite volume 
effects [IB] . The right panel in Fig. 1 shows the autocorrelation of the topological charge 
in the a = 0.12 and 0.09 fm ensembles. As expected, the autocorrelation time is larger 
for this quantity than for the others, and is much larger on the finer ensemble. The long 
autocorrelation time means that it will be important to check the size of finite volume effects; 
such runs are planned. 

In the sections below, the static quark potential was computed using block sizes of 50 
time units for a = 0.15 and 0.12 fm, and 60 time units for a = 0.09 fm. For the pseudoscalar 
meson plot, block sizes of 20 and 24 time units were used. Autocorrelations for the rho and 
nucleon masses are small, and were neglected here. 

IV. THE STATIC QUARK POTENTIAL 

Although it is not a physical observable, the potential between two infinitely heavy test 
quarks is well defined on the lattice and can be computed with high precision and com- 
paratively little effort. Therefore it has become conventional in lattice simulations to use 
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FIG. 2: The static quark potential with the HISQ and the asqtad actions. The HISQ results are 
from the a ~ 0.09 fm run, and the asqtad results are from a lattice with almost the same lattice 
spacing and light-quark mass about 0.2 times the correct strange quark mass (ami = 0.00465). In 
order to match the potentials, the plot is in units of n, while rulers in units of the lattice spacing 
are shown at r\V(r) = 0. A constant has been added to each potential so that V(r\) = 0. The 
solid lines (essentially superimposed) show the fit from Eq. ^ for the two runs, (evaluated with A 
set to zero). The inset magnifies a part of this plot at short distance to show the lattice artifacts 
discussed in the text. 
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a length scale based on the static quark potential to relate lattice simulations with differ- 
ent couplings, and to translate the dimensionless results of lattice simulations into physical 
units. We generally use r± defined by r\F{r\) = —1, where F{r) is the force — 9 v ^ . The 
scale r defined by TqF^q) = —1.65 is also commonly used. (The idea behind scales of this 
form [T7j is that they locate the transition region between the Coulomb potential at short 
distances, r 2 F(r) = — |a, and the linear potential at long distances, r 2 F(r) = —ar 2 .) 

In order to determine r\, we measure the static potential at discrete distances r 2 /a 2 = 
n 2 + n 2 + n 2 ,, and, for a range of r approximately centered at r±, we fit it to the functional 
form[18j 

B 



V(R) = C+- + aR + X 
K 











lot 





(3) 



Here j^\ lat is the the lattice Coulomb potential, ^\ lat = 47r J ^4sD^(p)e ipR , with D^lp) 
the free lattice gluon propagator calculated with the Symanzik improved gauge action, and 
1/R is the continuum Coulomb potential. 

Figure [2] shows the static quark potential at a w 0.09 fm for the HISQ ensemble and a 
corresponding asqtad ensemble. Overall, the two potentials are very similar. For reference, 
the value of r\ja for this HISQ ensemble in Table[l]came from a fit to the range a/5 < r/a < 6, 
or 0.63 < r/ri < 1.70. The inset in Fig. [2] makes visible some of the lattice artifacts at 
short distance. In particular, the HISQ point at rjr\ = 0.57 and the asqtad point at 0.53 
correspond to separation (2,0,0) along a lattice axis, and are visibly displaced below the 
trend. Note that artifacts of this kind are not decreased with the HISQ action, and we 
do not expect them to be decreased. In fact, in the continuum limit we expect them to 
be described by Eq. ^ with A = B. (The fit to this potential has B = —0.441(6) and 
A = —0.52(11).) Artifacts like this, at fixed number of lattice spacings, simply move to 
r = in the continuum limit. Also note that these artifacts diminish quickly with increasing 
r. For example, the HISQ point at rjr\ = 0.95 is really two points, with f/a = (3, 0, 0) and 
(2,2,1), and the difference between the two potential values is invisible. We expect that 
these short distance lattice artifacts in the static quark potential are mostly controlled by 
the gauge actions, which differ only in the fermion contributions to the one loop corrections. 

We do expect scaling violations proportional to a 4 and to a 2 a 2 at physical distances for 
both actions, and these would be visible in quantities like r Q /ri or r 1A /cr. However, it is not 
possible to make a definitive comparison of scaling violations in these quantities between the 
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TABLE III: Parameters of the potential fits in Fig. [2] As discussed in the text, in this comparison 
the fit ranges for the HISQ potential were chosen to match those used for the asqtad potential, and 
so these tabulated parameters differ slightly from those used in the rest of this paper. Note that 
the lattice mass is regularization dependent — in both of these ensembles the light quark mass is 
about one fifth of the correct strange quark mass. 





asqtad 


HISQ 


Fit range 


2.01-6.5 


2.01-6.5 


Time separations 


5-6 


5-6 


10/ 5 2 


7.085 


6.30 


ami/ am s /am c 


0.00465/0.031/na 


0.0074/0.037/0.440 


Ca 


0.849(3) 


0.824(3) 


B 


-0.432(4) 


-0.450(5) 


or? 


0.568(6) 


0.554(6) 


n/a 


3.697(7) 


3.510(7) 



two actions yet, since the addition of the dynamical charm quark to the HISQ simulations 
could also have small effects on these quantities. 

For reference, Table III shows the parameters of the fits in Fig. |2j defined in Eq. (j3|. 
Note that in this figure the fitting range used for the HISQ run is the same as used for the 
asqtad ensemble, and so differs from that used in finding the value of rq/a in Table [Tj Also 
note that the quantity ar\ parameterizes the potential in the range around rq, and should 
not be used as a measurement of the long distance string tension. Finally, note that since 
the dimensionful parameters are expressed in units of r±, which is found from the same fit, 
one relation between B and a is automatically enforced. In Fig. [2] this constraint forces 
both fits to have the same slope at r = rq (since rq is defined by the slope (force) at this 
distance), and a constant was subtracted to make both fits be zero at this point. 



V. SCALING TESTS 



Reduction of taste splittings among the pion masses with HISQ valence quarks was 
demonstrated with quenched gauge fields in Refs. [TJ |2], and with asqtad sea quarks in 
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FIG. 3: Taste splittings among the pions. The asqtad results used configurations with 2+1 flavors 
of dynamical quarks, and the HISQ results 2+1+1 flavors. The quantity plotted is r\ (M% — Mq), 
where M n is the mass of the non-Goldstone pion and Mq is the mass of the Goldstone pion. 
Reading from top to bottom, the non-Goldstone pions are the tt s (box), ttq (fancy box), 7Tj (fancy 
plus), iTio (plus), 7Ty (diamond), 7Tj5 (cross) and 7ros (octagon). r\ (M% — M G ) is known to be 
almost independent of the light-quark mass. The vertical bar at the upper left shows the size of 
a factor of three, roughly the observed reduction in taste splittings, while the sloping solid line 
shows the theoretically expected dependence on lattice spacing. Nearly degenerate points have 
been shifted horizontally to improve their visibility. 
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Ref. [3], and there is little reason to expect it to be different with dynamical HISQ sea 
quarks. However, in view of the importance of this quantity, we show splittings for all of the 
different tastes of pions in Fig. |3j comparing results with HISQ quarks (both valence and 
sea) to earlier results with asqtad quarks. In this figure, we see that the expected reduc- 
tion in taste splittings happens, with roughly a factor of three reduction relative to asqtad 
calculations at the same lattice spacing. 

The main purpose of this study was to see if the improvements in the action designed 
to reduce taste symmetry violations translate into decreased lattice spacing dependence in 
other quantities. We begin with the mass of the light-quark vector meson, or p. In Fig. [4] 
we show the mass of the p meson in units of r\. Here we have asqtad results for several 
light-quark masses at each lattice spacing, but HISQ results for only one light-quark mass. 
The light-quark masses themselves are regularization dependent, so to plot asqtad and HISQ 
results on the same footing we use the Goldstone pion mass in units of r\ for the horizontal 
axis. Note that for m\ = 0.2m s , the light-quark mass used in the HISQ simulations, and 
for the lattice sizes used here (< 2.9 fm), the vector meson is stable against decay to two 
pions. Results for the nucleon mass are similar, and are shown in Fig. [5] In Figs. [4] and [5] 
the HISQ masses show smaller dependence on the lattice spacing than the asqtad masses, 
with the same continuum limits within the statistical errors. Roughly speaking, the HISQ 
results are similar to the asqtad results at the next smaller lattice spacing. 

Although we have chosen to present these results as improved scaling of the p and nucleon 
masses, since we are plotting the dimensionless quantities M p r 1 and M N r 1 , they could equally 
well be described as improved scaling of r\ when a hadron mass is chosen to be the length 
standard. 

The pseudoscalar meson decay constants are important for lattice determinations of CKM 
matrix elements, and can be computed with high precision. In fact, our current best de- 
termination of j*i in physical units comes from matching the asqtad lattice results to the 
physical value of /„.. These decay constants for light quarks have been extensively studied 
using the asqtad ensembles [U [20]. The HPQCD collaboration has computed these decay 
constants in a mixed action calculation, with HISQ valence quarks on the asqtad sea quark 
ensembles, and used them in a determination of the physical value of r\ |21j . Figure [6] shows 
the pseudoscalar decay constant with one of the valence quarks fixed at approximately the 
strange quark mass as a function of the mass of the other valence quark. (At the physical 
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FIG. 4: Vector meson {p) masses in units of r\. Here the bold (red) points are the HISQ simulations 
with mi = 0.2m s , and the lighter (blue) points are asqtad results for various light quark masses. 
The a ~ 0.06 fm asqtad point immediately to the right of the a ~ 0.09 fm HISQ point has been 
displaced to the right to make it visible. It in fact falls on top of the a ~ 0.09 fm HISQ point. 
The cross sign at lower left is the physical p mass. The error on the physical mass point is just the 
error on the physical value of r\. 
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FIG. 5: Nucleoli masses in units of r\. Here the bold (red) points are the HISQ simulations with 
m/ = 0.2m s , and the lighter (blue) points are asqtad results for various light quark masses. The 
cross at lower left is the physical nucleon mass. The solid magenta line is a continuum extrapolation 
of a chiral perturbation theory fit to the asqtad nucleon masses, while the dotted green lines are 
from the same fit at finite lattice spacing [19J. 

light-quark mass, this is just fx-) To facilitate the comparison, we have used the ratio of 
the light quark mass to the corrected strange quark mass in the corresponding ensemble for 
the horizontal axis. The reduction in lattice artifacts is obvious, and it can also be seen that 
the HISQ points lie near the continuum limit of the asqtad points. Once again, we remark 
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FIG. 6: Pseudoscalar decay constant. One valence-quark mass, tua, is varied while the second is 
held fixed near the strange-quark mass. All ensembles used a light sea quark mass of about 0.2 
times the strange-quark mass. The left hand panel shows asqtad results for four different lattice 
spacings and the right hand panel shows HISQ results for three lattice spacings. 



that since the plotted quantity is rifps, this could equally well be described as scaling of r\ 
or scaling of fps- 

The topological susceptibility is a particularly important test here, since it is computed 
solely from the gluon configurations that are generated, without involving HISQ or asqtad 
valence quarks. Therefore improvements in the scaling of the topological susceptibility di- 
rectly test whether the change of the sea-quark action has the expected effect on the gluon 
configurations that are generated. Our technique for calculating the topological suscepti- 
bility is set out in detail in Ref. [23 J. Here we just note that this technique is based on 
measurement of a density-density correlator, and hence is not limited by long autocorrela- 
tion times for the overall topological charge. Figure [7] shows the topological susceptibility 
for most of the asqtad ensembles, and HISQ results for the a ~ 0.12 fm and a ~ 0.09 fm 
ensembles. The HISQ point with a 0.12 fm lies near the asqtad points with a ~ 0.09 
fm, and the HISQ point with a « 0.09 fm is near the asqtad points with a ~ 0.06 fm, 
demonstrating a decrease in lattice artifacts. Note that the HISQ points are to the left of 
the corresponding asqtad points, which are indicated by arrows in the figure. This is because 
the horizontal axis is the mass of the taste singlet pion (the heaviest pion taste), and the 
reduction in taste symmetry breaking moves the points to the left. It is the movement down 
relative to the asqtad points that represents an improvement in the gluon configurations. 
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TABLE IV: Lattice spacings in fm from r\ = 0.3117 fm, f ss with asqtad valence quarks, and f ss 
with HISQ valence quarks. The first five columns identify the ensemble by the sea-quark action, 
the gauge coupling 10 /g 2 , and the sea-quark masses. The horizontal line separates ensembles with 
asqtad sea quarks (above) from those with HISQ sea quarks (below). The values for HISQ valence 
quarks on asqtad sea ensembles are taken from Ref. |21| . The errors on a(r\) are statistical only — 
they do not include the errors in n = 0.3117(6) (tH) fm - Similarly, the errors on a(f ss — asqtad) 
and a(f ss — HISQ) for the HISQ ensembles do not include any errors in the physical value of f ss . 



The numbers following the f ss lattice spacings are the value of the valence strange quark mass am s 
at which the desired ratio is obtained. We use the values f ss = 181.5 MeV and f ss /M ss = 0.2647 
from Ref. [2]. 



Action 10/ g 2 ami am s am c 


a(ri) 


a(fss — asqtad) am s 


a(f ss - HISQ) am s 


asqtad 6.76 0.01 0.05 
asqtad 7.09 0.0062 0.031 
asqtad 7.46 0.0036 0.018 
asqtad 7.81 0.0028 0.014 


0.1178(2) 
0.0845(1) 
0.0588(2) 
0.0436(2) 


0.1373(2) 0.0467 
0.0905(3) 0.0286 
0.0607(1) 0.0187 
0.0444(1) 0.0133 


0.1264(11) 0.0553 
0.0878(7) 0.0362 
0.0601(5) 0.0233 
0.0443(4) 0.0163 


HISQ 5.80 0.013 0.065 0.838 
HISQ 6.00 0.0102 0.0509 0.635 
HISQ 6.30 0.0074 0.037 0.440 


0.1527(7) 
0.1211(2) 
0.0884(2) 


na na 
na na 
na na 


0.1558(3) 0.0720 
0.1244(2) 0.0549 
0.0900(1) 0.0374 



VI. USING f ss TO SET THE SCALE 



In Figs. |4]-[6]it can be seen that rifx and the hadron masses in units of r\ all increase as 
the lattice becomes coarser. This common dependence on lattice spacing could be absorbed 
into a lattice spacing dependence of r\. Put more simply, we could use one of these quantities 
to set the lattice spacing. Such a procedure has been introduced and studied by the HPQCD 
collaboration in Ref. [21]. In particular, they use the decay constant of a fictitious "unmixed 
ss" pseudoscalar meson, which is an isospin non-singlet meson with both valence quarks 
having mass m s , to set the scale. We call this decay constant f ss . Like r\, f ss is not a 
quantity that can be directly determined from experiment, and so, like r\, its physical value 
is eventually determined by matching to some precisely known quantity such as f n or mass 
splittings of heavy quark mesons. In practice, the HPQCD collaboration determines f ss and 
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FIG. 7: The topological susceptibility. Points with the asqtad action are shown for several 
lattice spacings and quark masses, and the HISQ results for a ~ 0.12 fm and a ~ 0.09 fm with 
mi = 0.2m s . For the horizontal axis we use the mass of the taste singlet pion, since in lowest order 
chiral perturbation theory the topological susceptibility is a function of this mass [22J . The curves 
in the figure come from a chiral perturbation theory fit to the asqtad data. The asqtad results 
are updated from Ref. [21] and are discussed further in Refs. [H [23]. The two arrows indicate the 
locations of asqtad points with lattice spacing and quark mass similar to the two HISQ points. (In 
the case of the a ~ 0.09 fm HISQ point, the quark mass falls between two of the masses of the 
asqtad points.) 
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the corresponding meson mass M ss using a next-to-leading-order chiral fit (augmented with 
discretization corrections) to their lattice data, and inputs of the experimental values for f„, 
fx, M w and Mk- In fact, lowest order chiral perturbation theory and these experimental 
values alone (without lattice data) gets within ~ 1% of the HPQCD results [21J. We prefer 
not to input the experimental value fx in such determinations, since we take fx as an 
output of our lattice calculations that gives a result for V us [UEO]- Indeed, f n , M w and M K 
alone are adequate for determining the physical scale and the quark masses mi and m s , and 
hence all light-quark quantities. 

An advantage of using f ss to set the scale on a given lattice ensemble is that it can 
be determined to high accuracy in the simulations. However, it has the disadvantage that 
it depends on the choice of valence quarks, so the lattice spacing assigned to a particular 
ensemble will depend slightly on whether it is determined with asqtad quarks, HISQ quarks, 
or some other formalism. 



Table [TV] shows the lattice spacings of the three HISQ ensembles used in this paper, and 
some comparable asqtad ensembles, using r± and / ss as the length standards. For the asqtad 
ensembles, we show the effect of using either asqtad or HISQ valence quarks to determine 
the lattice spacing. Note that, as expected, the differences among the scale determinations 
decrease as the lattice spacing decreases. The value of / ss and the corresponding strange 
quark mass am s were determined by fitting a quadratic polynomial through masses and 
decay constants at valence masses equal to 1.0, 0.8 and 0.6 times the sea strange quark mass. 



Table [TV] also shows the value of the strange quark mass am s given by this interpolation or 
extrapolation. Figure [8] shows the differences in length scale (relative to the determination 
from as a function of lattice spacing. In this figure it can be seen that these differences 
are vanishing in the expected way as the lattice spacing decreases. 

In Fig. [9] we show the rho mass data from Fig. [4] replotted using / ss to set the scale. 
Replotting the nucleon masses in Fig. [5] would produce similar results. (Of course, one could 
then make a plot showing the dependence of rif ss on sea-quark mass and lattice spacing.) 

VII. CONCLUSIONS 

Using simulations with a fixed and unphysically large light-quark mass, we see that di- 
mensionless ratios of several hadronic quantities show smaller dependence on lattice spacing 
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FIG. 8: Differences in determinations of the length scale using different standards. In the legend, 
the symbol types are labelled as "valence on sea". The "HISQ on asqtad" points are taken from 
Ref. [21]. 

with the HISQ action than with the asqtad action. Roughly, for the quantities that we 
checked, HISQ simulations at lattice spacing a appear to have similar lattice artifacts as 
asqtad simulations at lattice spacing |a, leading to substantial savings in simulation costs. 
This program is continuing with computations at different light sea-quark masses, so that 
both the extrapolation to the continuum limit and the extrapolation to the physical light- 
quark mass can be controlled. 
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FIG. 9: Vector meson (/?) masses in units of f ss . The data and the meaning of the symbols are the 
same as in Fig. |4j The vertical and horizontal scales in the figure correspond to the same ranges 
as in Fig. |4j 
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Appendix A: Gauge and fermion actions 

For completeness, we summarize the gauge and fermion actions in this appendix. 

The gauge action is a tadpole-improved [25] one-loop-Symanzik-improved gauge action 
[26] including the effects of the quark loops in the one-loop coefficients [8]. The number of 
flavors rif is set to four in these simulations. 

This gauge action involves three kinds of loops: the lxl loop, or plaquette P, 
the 2x1 loop, or rectangle R, and the twisted loop T, which traverses paths such as 
+x, +y, +z, —x, —y, —z. Then 

Sg = P (CP E ( X " ^ Re Tr ( P )) + C « E i 1 ~ l Re Tr ^)) + Ct E i 1 ~ l Re Tr ( T )) ) ' 

(Al) 

where the sums run over all distinct positions and orientations of the loops. The coefficients 
are 

C P = 1.0 

°R = -(0-6264- 1.174671/) ln(wo)) 

C t = \ (0.0433 - 0.0156ra/) ln(u ) . (A2) 
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In this expression the strong coupling constant appears in the form a s = — ln(ito)/l -303615. 
With this normalization (3 = p. We determine the tadpole coefficient uq from the average 
plaquette, uq = ((ReTrP) /3) 1 / 4 . 

The fermion factor in the partition function is 

\n(S f ) = H(det(2p + 2m f )f 4 . (A3) 
/ 

The Dirac operator p is constructed from smeared links. Two levels of smearing are used, 
with a projection onto an element of U(3) after the first smearing. The fundamental gauge 
links are U^x), the fat links after a level one fat7 smearing are V^(x), the reunitarized 
links are W^x), and the fat links after level two asqtad smearing are X^{x). The first level 
smeared links V are constructed from the U as a sum over products of links along paths 
from x to x + fi, or parallel transports. 

W = 52Y[uM(path) (A4) 

paths path 

Table [V] gives the coefficients used in the two levels of smearing. The nearest neighbor 
part of p uses the twice-smeared links X while the third nearest neighbor part uses the 
once-smeared and unitarized links W: 



2p x>y = Y^ { <W,2/ X M ( x ) - S x- ^y X t (x-fi)} 
+ (1 + e N ) {5 x+ ^ y W^{x)W^x + ftW^x + 2/2) 

- 5 x . 3(l , y Wl(x - 3fi)Wt(x - 2jl)Wl(x - A)} • (A5) 



In Eq. (|A§ and Table |V| 6n is a mass-dependent correction to the tree-level improvement 

of the quark dispersion relation, or the "Naik term." This correction is negligible and set to 

zero for the light and strange quarks. For the charm quark we use 

27, , 2 327 , , 4 15607 , , B 73697 . ,„ . 

ejv = am, H am, am, am, . (A6) 

40 v ; 1120 v ; 268800 v 1 3942400 v 1 y ' 

In this expression am c is the bare mass in the quark action, and this formula combines 
Eqs. (24) and (26) in Ref. [3]. The numerical values of ejv used in our simulations are given 
in Table [H 
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TABLE V: Paths and coefficients used in smearing the links. It is understood that all distinct 
rotations and reflections of each path are used in the action. In specifying paths in this table, 
directions x and y, etc. are different. The multiplicity is the number of such paths contributing to 
a single smeared link. The first block of the table gives the coefficients used in the "fat 7" smearing 
used to construct V from U. The second block gives the coefficients in the "asqtad+" smearing 
used to compute X from the unitarized links W, and the final line is the coefficient of the third 
nearest neighbor term. Note that the coefficient of the "Lepage" term that corrects the form factor 
at small momenta is twice that of the single smearing asqtad action. 



Name Path 




Multiplicity Coefficient 


Single link +x 




1 


1/8 


3-staple +ij + x — y 




6 


1/16 


5-staple +y + z + x — z 


- V 


24 


1/64 


7-staple +y + z + t + x 


-i - z 


-y 48 


1/384 


Single link +x 




1 


1/8 + 3/4 +1/8(1 + cat) 


3-staple +y + x — y 




6 


1/16 


5-staple +y + z + x — z 


- V 


24 


1/64 


7-staple +y + z + t + x 


-t - z 


-y 48 


1/384 


"Lepage" +y + y + x - y 


- V 


6 


-1/8 


"Naik" +x + x + x 




1 


-1/24(1 + 6*) 



Appendix B: HISQ force calculation details 

Here we summarize the details of the HISQ force calculation. For clarity in this and the 
next appendices we suppress the x-dependence and the direction index in the notation of 
the links. Most of this material appeared earlier in Ref. [ID] . For the force calculation we 
adopted the strategy of Refs. [27] and [28] , in which the derivative of the smeared action is 
calculated by repetitive application of the chain rule: 

dS f _ dS f OX dW dV 

W ' ~ ~dX ~dW ~dV dU ' ( ' 

where Sf is the fermion part of the action, U are fundamental gauge links, V, the fat links 
after level one fat 7 smearing, W, the reunitarized links, and X, the fat links after level 
two asqtad smearing. In our code, for the parts that involve smearing, we follow the same 
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FIG. 10: Time history of the maximum (over lattice sites) magnitude of the fermion force and the 
minimum determinant of the fattened links after the first level of smearing (V). This exploratory 
run was done on a 20 3 x 64 lattice at j3 = 6.75, ami = 0.2am s , am s = 0.05, am c = 0.6 and uq = 0.9. 
This approximately corresponds to the a ~ 0.12 fm ensemble in Table [I] at f3 = 6.0. The difference 
in f3 is due to the use of a different gauge action in the earlier studies. 



procedure as for the asqtad action. This procedure is described in Refs. [29] and [3"0] . so we 
do not repeat it here. The algorithm for the reunitarization part is detailed below. 

We have chosen to project links to U(3), rather than SU(3) as in the original 
HPQCD/UKQCD formulation, for two reasons: 
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1. SU(3) projection requires calculation of the third root of the determinant, which in- 
volves a phase that can initially be restricted to, e.g., the interval [— 7r/3, 7r/3). How- 
ever, during the molecular dynamics evolution, this phase has to evolve continuously 
(to prevent the appearance of 5 function-like forces) and may cross into [7r/3,27r/3) 
interval, and so on. Thus, SU(3) projection requires tracking the evolution of the 
phase for each link during molecular dynamics. 

2. For the U(3) group, different methods of projection yield the same answer for the 
projected link, W. 

For instance, the default method in our code is polar projection: one builds a Hermitian 
matrix 

Q = 1/V (B2) 

and then 

W = VQ- 1/2 (B3) 

belongs to U(3), i.e., 

W^W = (Q- l/2 ) ] V ] VQ~ l/2 = Q- l l 2 QQ' 1 / 2 = 1. (B4) 

It is important that closed form expressions for Q~ x l 2 can be derived [31J and, thus, the 
whole procedure can be implemented analytically. 

One may expect, given that Q" 1 ^ 2 is a singular operation, that when one of the eigenvalues 
of Q is close to 0, the numerical accuracy in evaluation of W becomes poor. In fact, in 
simulations one occasionally encounters large deviations from unitarity: 

\W ] W 0{1). (B5) 

Such situations are rare, but the contribution of (systematic) errors of this kind can be large. 
Therefore, we also implemented the singular value decomposition (SVD) algorithm, which 
is slower, but is used only in exceptional cases. We decompose: 

V = AEB\ (B6) 

A,Be U(3) and E is a positive, diagonal matrix. (The values on the diagonal are called 
the singular values of V .) Then we have, simply, 

W = AB ] . (B7) 
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It is easy to see that (B3) gives the same result as (B7): 



Q = V ] V = BTjA^ATjB^ = BY,EB ] = BEB ] BEB ] = (BEB ] ) 2 , 
Q- 1 ' 2 = (fiSfit)- 1 = (fit)- 1 s-i^-i = ss-ifit ; 

W = AY,B ] BY 1 - 1 B ] = AB ] , 

as it should be. The SVD algorithm of Golub and Reinsch [32] is numerically stable even in 
the case of exactly zero eigenvalues. 

Another popular projection algorithm is "trace maximization" [33]: find W G U(3) such 
that it maximizes 

VL = ReTi{V ] W) . (B8) 

Let us again use SVD on V: 

V = AEB ] Q = ReTr {BEA^W} = ReTr {EA^WB} . (B9) 

Since E is positive, clearly Q is maximized when 

A^WB = 1 W = AB ] (BIO) 



and we arrive at (B7) again. In the SU(3) case an extra phase present in (B9) would lead to 



W different from (B7). To summarize, we use the polar projection (B3) replaced by SVD if 
small eigenvalues of Q are encountered. 

Appendix |c| gives the details of calculation of Q~ 1 ^ 2 and its derivative dQ^ 1 ^ 2 / dV within 
the approach of Refs. |TTj and [3T] based on the Cayley-Hamilton theorem. 

Let us now turn to the calculation of the force. During the molecular dynamics evolution, 
one encounters (more often on coarser ensembles) matrices V that have small eigenvalues. 
Let us consider the U(l) group for simplicity. Then V is just an arbitrary complex number 
V = re 10 . The projection onto U(l) is W = e td . The derivative that enters the force 
calculation is 

dW _ fdW\ d(W,V^) d(W,V^) d(r,6) 1 . 

(Bll) 



dV \dV J vi 9(V,Vt) d(r,6) d(V,V^) 2r 

Thus the derivative is inversely proportional to the magnitude of V, or, in the U(3) case, 
to the smallest singular value of V (or eigenvalue of Q), which is not protected from being 
zero. Thus, on rare occasions one has to deal with exceptionally large forces that give large 
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contributions to the action, but originate from a single link. In Fig. 10 we show the evolution 
of the minimal det \ V\ over the lattice and maximal value of the norm of the fermion force 
on a logarithmic scale. One can easily see the correlation: the lower det \ V\, the higher the 
force. 

To circumvent this problem we introduce a "cutoff" in the force calculation by replacing: 



W = VQ 



-1/2 



W = V(Q + SI) 



-1/2 



(B12) 



where / is the unit matrix, whenever the smallest eigenvalue of Q is less than 5. In the 
ensembles used in this paper we set 5 = 5 x 10~ 5 . In tuning 5, we weigh two competing 
issues: The value of 5 should be large enough to suppress an exceptional contribution from 
a link, but small enough not to modify too many forces on the lattice. If 5 is too large, the 
evolution will be smooth, but the fluctuation of the action will be large, usually leading to 
rejection of the trajectory. Note that we modify W only in the force calculation, and we 



use the original Eq. (B3), or Eq. (B7) for nearly singular matrices, to calculate the action at 



the accept/reject step. That is, the modification (B12) amounts to using a different guiding 



Hamiltonian during the evolution, while the Metropolis step insures the desired distribution. 



Appendix C: Algebra for reunitarized links and their derivatives 



To make this presentation self-contained we include Eqs. (Cl)-(Cll) from Hasenfratz, 



Hoffmann and Schaefer [31], preserving the notation of the original. 



The inverse square root of a non-singular matrix Q entering Eq. (B3) is given by the 
Cayley-Hamilton theorem as a polynomial of Q: 



Q 



-1/2 



fo + fiQ + f 2 Q 2 



(Ci) 



Since Q is Hermitian, it has nonnegative eigenvalues that can be found by solving the 
characteristic equation 



where 



g 3 - c g 2 - ( ci - -eg j g - (c 2 - c ci + ^ ) = 0. 



„ ^^^-trQ" +1 , n = 0,l,2. 



n + 1 



(C2) 



(C3) 
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The solution of the cubic equation (C2) is 

c 



9r, 



3 + {n ~ 1) 



2vr\ 
"3"/ 



I , n = 0,l,2, 



where 



(C4) 



(C5) 



It is convenient to define the symmetric polynomials of the square roots of the eigenvalues 



w 



y/9o9i92 ■ 



(C6) 



In the diagonalized form the expression (CI) can be rewritten as an equation for /j 

1 1 '/(i '/(! ./n '/(, 

-1/2 



which has the solution 



1 go go 
1 9i g\ 

\ i #2 0! y 



./o 

/l 

V /2 / 



01 

-1/2 
V^2 J 



fo 

h 

h 



—w{u 2 + v) + UV 2 

w{uv — w) 

—w — u 3 + 2uv 

w(uv — w) 
u 



w(uv — w) 



The derivative dfi/dcj can be written as 



dfi dg k 



13 dc 3 d 9k d Cj 



After rescaling ( C9 ) by the common denominator 



(C7) 



(C8) 



(C9) 



Qj = dB i:j , d = 2w 3 (uv -w f 



(CIO) 
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a closed-form expression for the symmetric matrix has been derived in Ref. [3T] 
C 00 = -w 3 u e + 3vw 3 u 4 + 3v 4 wu 4 - v 6 u 3 - Aw 4 u 3 - I2v 3 w 2 u 3 

O Q Q ^ O A A O ^\ 

+16f w u + 3u m — Svw u — 3v w u + w + v w , 
Cqi = -w 2 u 7 - v 2 wu 6 + v 4 u 5 + 6vw 2 u 5 - 5w 3 u 4 - v 3 wu 4 - 2v 5 u 3 

Q Q q Q O A Q /I Q O O Q 

— 6f tt) u + 10f ty « + 6f «m — 3u; it — 6n to u + 2v w , 

O ^ O/l /I Q O Q Q O Q O 00 Q 

602 = w u + v wu — v u — Avw u + Aw u + 3v wu — 3v w u + m , 
C u = -wu 8 - v 2 u 7 + 7wu7u 6 + Av 3 u 5 - 5w 2 u 5 - lQv 2 wu 4 - Av 4 u 3 + \Qvw 2 u 3 

-3w 3 u 2 + I2v 3 wu 2 - Y2v 2 w 2 u + 3vw 3 , 
Ci 2 = wu 6 + v 2 u 5 - 5vwu 4 - 2v 3 u 3 + Aw 2 u 3 + 6v 2 wu 2 - 6vw 2 u + w 3 , 
C 22 = -wu 4 - v V + 3vwu 2 - 3w 2 u . (Cll) 

In the following, differentiation with respect to V at fixed V' is always assumed. We use 
explicit color indices to show how different contractions and direct products of matrices are 
built. 

The derivatives that enter in the calculation of the fermion force are 

m; ~ 9v kl - 5M >« + Vim w M ' (C12) 

awl WQ- 1/3 )*JO _ d(Q- 1/2 W T/t 

v mj . \S j1 - 6 ) 



Also, 



dV kl dV kl dV kl 

dQ lj 



dV kl 

The central component of the calculation is 

0(Q- 1/a )« . d 



VlSir (C14) 



{foSij + hQij + f2(Q 2 )ij) 



9 fo r , dfl „ , r r c , ^ ,^ 2 N 



) (C15) 



9Qpq 9Qpq 9Qpq 



From the definition (C3) it follows that dc n /dQ pq = (Q n ) pq . Then 
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We define 



dfo ■ 



R 



qp 



Sqp — 



dfl 

9Qpq 

d_h 

9Qpq 



qp j 



BloS qp + BiiQ qp + Bi 2 (Q z 

■ B 2 o5 qp + B 2 iQq P + B 2 2(Q 2 ) 



iqp , 



qp ■ 



(C17) 
(C18) 
(C19) 



Substituting (C17), (C18) and (C19) into (C15) and Eq. (C15) in (C12) and (C13) we obtain 



the final result 



dVi 



ki 



= S^Q- 1 ^ + [fi(VV% + f 2 (VQV%] <J y + f 2 (VV%Q lj 

+ V^PV^u + (VQ) tj (RV% + (VQ%(SV% , 

= 'fiVl + f 2 (QV%] V* + f 2 Vl(QV% 

+ Vt(PV% + (QV%(RV% + (Q 2 V%(SV^) 



(C20) 



ik 



(C21) 



The calculation of the fermion force from the reunitarized links proceeds as follows: 



1. The eigenvalues of the Hermitian matrix Q are calculated with Eq. (C4). 



2. det \Q\ is compared with the product gog\g 2 - If the relative error is larger than 10 -8 or 
any eigenvalue is smaller than 10 -8 the singular value decomposition of V is performed 
and the eigenvalues are set to 



gi = a i , i = 0, 1, 2, 



(C22) 



where Oi are the diagonal elements of the matrix E in Eq. (B6). 



3. Additionally, if any of the eigenvalues is smaller than (an adjustable parameter) 5 
5 x 1CT 5 the eigenvalues are modified to 



9i ->■ 9i + 8. 



(C23) 



(This corresponds to the force "cutoff" in Eq. (B12). 



4. With these eigenvalues the coefficients fi and the elements Bij are calculated from 



Eq. (C8) and (Cll). 



5. Finally, the force is calculated from Eq. (C20) and (C21) 
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In the MILC code we have also implemented two other methods for calculating Q 1//2 
and its derivative: 



1. a rational function approximation, 

2. an iterative evaluation of Q^ 1 ^ 2 with the derivative replaced by finite difference. 



We found that the analytic evaluation via Eq. (C20) and (C21) is superior to the other 
methods due to its higher precision and speed. 

Appendix D: Treatment of the HISQ charm quark 

The tree-level discretization errors are 0((ap M ) 4 ) and are negligible for light quarks. 
However, at the lattice spacings listed in Table |TJ the charm quark mass is in the range 
am c ~ 0.4 — 0.8 and therefore the discretization errors are larger. The leading tree-level 
0((am c ) 4 ) error can be removed by retuning the coefficient of the third-nearest-neighbor 
(Naik) term [3], using the expansion in Eq. (A6). As can be seen from Table [TJ at the finest 



lattice, a m 0.09 fm, it is quite small, = —0.120471. The effect of the correction in 



Eq. (A6) has been studied in Ref. [3]. To check the quality of charm quark physics in our 



ensembles, we computed the speed of light for the 1] C meson by calculating its propagator at 



several non-zero momenta. The result is shown in Fig. 11, where 



and the momenta are rescaled by the lattice size L s 

n 2 =p 2(Ls\\ (m) 



2n / 

For the finest a ~ 0.09 fm ensemble the error in the dispersion relation is below 2%. As 
expected, these r] c dispersion relations are very similar to those found for HISQ valence 
quarks on asqtad sea quarks, shown in Table V of Ref. [3]. 

For dynamical simulations, the mass-dependent correction to the Naik term requires 
the use of different sets of smeared links for light quarks and the charm quark. Since the 
difference in the Naik term enters at the second level of smearing, it is advantageous to 
regroup the force calculation as described in the following. Let denote the fat links 
after level two asqtad smearing for the light quarks, for which is set to zero, and 
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FIG. 11: The speed of light for the r/ c meson calculated at several values of momenta. 



denote the fat links for the charm quark. Then the fat links for the charm quark can be 
written as 

X (c) = x (o) + £nAX ; ( D3 ) 

where AX contains only one-link and three-link paths. (This can be easily seen from Ta- 
ble [Vj) The derivative is 



ox^ ax^ oax 



dW " dW 



+ e N - 



dW 



(D4) 



The fermion force in Eq. (Bl ) contains contributions from the light (u : d and s) quarks, and 
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from the charm quark. 

dS f dS f dXW dW dV dS f dX® dW dV 

~dU ~ dX(°) dW ~dVdU + dX^ c ) dW ~dV dU ' ^ ' 

The calculation of the force for multiply smeared actions proceeds from the last level of 
smearing to the first one. Therefore, operations with X^ and X^ links are done first and 
can be combined before the reunitarization part: 

dS f f dS f dX^ dS f dX^\ dW dV 
~dU ~ \dX(°) dW + dX( c ) dW ) ~W dU 

f dSf dS f \ dX® dS f dAX\ dW dV 

\dx^ + dx^) ) ~dw~ + eN dx^) ~dw~ ) WdU' ( ' 

After the dAX/dW contribution is separated, the number of operations needed for the 
HISQ fermion force is reduced to slightly more than twice the number needed for the asqtad 
fermion force. This is because the most time-consuming part of the calculation is related to 



3-, 5- and 7-staple paths that have high multiplicity. In our final form (D6) they are present 
only in dX^/dW and dV/dU. 
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